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Abstract. This work develops new numerical methods for the solution of the tomography 
problem in domains with reflecting obstacles. We compare the solution's performance for 
Lambertian reflection, for classical tomography with ubroken rays and for specular reflection. 
Our numerical method using Lambertian reflection improves the solution's accuracy by an 
order of magnitude compared to classical tomography with ubroken rays and for tomography 
in the presence of a specularly reflecting obstacle the numerical method improves the solution's 
accuracy approximately by a factor of three times. We present efficient new algorithms for the 
solution's software implementation and analyze the solution's performance and effectiveness. 



1. Introduction 

Let f(x) be a continuous function in f2, where Q = fio\^i an d is a compact convex set 
in R 2 with a smooth boundary and Ox a convex obstacle with a smooth boundary such that 
S7i C Oo C K 2 . Consider dQo as the observation boundary of f^o and the points of this boundary 
as transmitters and receivers of ray signals or ray solutions of the wave equation 

u tt -c 2 (x)Au = (1) 

where c(x) > is the variable speed of sound in Q and u\qq 1 = 0. 

The data for the tomography problem are all integrals J f{l)dl = C 7 where 7 are rays in 
Vi that have both of their end points in 8£Iq. One end point is a transmitter and the other 
end point of each of the rays a receiver. The classical tomography problem is to find f(x) in £1 
knowing the values of all integrals C 7 where 7 are all straight line segments or unbroken rays 
in Q. In the case of a domain Qq without an obstacle C this problem is widely studied 
theoretically and numerically [1, 2, 3, 4]. When there are obstacles present, this problem is 
much less studied. Key theoretical work is done in [5, 6, 7, 8] for domains with one obstacle and 
with both broken rays, i.e. rays reflecting at the obstacle, and unbroken rays. A key result from 
[5] is that the tomography problem in the presence of one reflecting obstacle is well-posed. In 
other words, if an obstacle C is present then we have the well-posed problem of recovering 
f(x) in ilo\^i from the set C 7 where 7 are all broken and unbroken rays in the domain starting 
and ending at the observation boundary. This problem is called the Broken Ray Tomography 
Problem [5]. 



In a basic tomography setup transmitters and receivers of wave signals are placed at the 
domain's boundary <9f2o- Ray signals are generated by the transmitters and received by the 
receivers. Travel times for signal propagation from transmitters to receivers are measured and 
these travel times T(A, B) are the values of line integrals of a function f(x) = where c(x) > 
is the speed of sound at point x € fio\^l- This measurement procedure gives the C 7 data for 
solving the Tomography Problem by relating signal travel times to the values of line integrals of 
f(x). Given sufficient data, f(x) and from here the velocity c(x) are computed with tomographic 
reconstruction algorithms [1, 2, 3, 4]. The tomographic algorithms presented in the next section 
require and compute the ray path 7 of each ray and, in order to compute these ray paths, we 
consider reflection models at d£l\ and models of the speed of sound in Q. 

The theory of broken ray tomography from [5, 6, 7, 8] is based on a reflection model at 9f2i 
that is mirror-like i.e. the angle of incidence is equal to the angle of reflection. In this paper, we 
consider a Lambertian reflection model in which incident rays are reflected at the obstacle in all 
possible directions and present results that show that the broken ray tomography reconstruction 
error is smaller when we consider Lambertian reflection. 

We consider a mathematical model of the speed of sound 

c{x) = c + e(x) (2) 

where x £ £1 that models the speed of sound as a continuous function close to a known constant 
c Q . It is shown in [9] that for sufficiently small e(x) waves propagate along the known geodesies 
of c D (x) when 

c(x) = c (x) + e(x) (3) 

where c(x), c (x) and e(x) are continuous functions in fi. The acoustic geodesies for constant 
speed of sound c (x) = c Q are straight lines when there is no obstacle. Therefore, for the model 
2, we consider two cases. In the first case, 7 = 71 is an unbroken ray composed of a straight 
line segment 71 . In the second case of a broken ray, 7 = 71 (J 72 is the union of two straight line 
segments that intersect at a reflection point at the obstacle. 
For unbroken rays the travel time or time of flight is 



J 7 c + e(x(s)) J fl " 



c + e(x(s)) 

and this model leads to the classical tomography problem with f(x) - 



c +e(x) ' 



For broken rays and a known obstacle, we know that the acoustic wave u(x) propagates along 
the known straight line segments 71 and 72. 71 and 72 are known because in addition to the 
time of flight, our data measurement procedure gives the end points of the ray 7 and its initial 
velocity, which in turn imply the reflection point of 7 at the known obstacle Q±. Then 

r^, . „s f ds f ds f ds 

+ e{x(s)) y 7 - c + e{x{s)) J^ 2 c 



+ e(x(s)) 

and, when data of this type is added to the set of measurements for the time of flight for unbroken 
rays, this gives the set C 7 for the Broken Ray Tomography Problem with f(x) = - _^ € ^ ■ In 
other words, the data set C 7 for the Broken Ray Tomography Problem contains the travel times 
of all broken and unbroken rays in the domain that start and end at the observation boundary. 



2. Numerical Solution of the Broken Ray Tomography Problem for Lambertian 
Reflection 

The first algorithm for the Broken Ray Tomography Problem is presented in [10] for a known 
obstacle Qi, specular reflection and the model 2 for the speed of sound. This work extends 
the first numerical solution of the Broken Ray Tomography Problem and develops an algorithm 
for finding the velocity structure of £1 for Lambertian reflectance at dQ±. The broken ray 
tomography problem for an obstacle with Lambertian reflectance is solved with the following 
algorithm that constructs the finite set of broken and unbroken rays and computes the associated 
ray travel times by numerical integration or from ray travel time data. 
Require: Domain 
Require: Obstacle f^i C Oo 

Require: Finite set of receiver points R = (x r ,y r ) on 80,q 

Require: Finite set of transmitter points T = (xt,yt) on d£lo 

Require: Finite set of obstacle boundary points H = (xh,Uh) on 

Require: Number of broken rays n&. 

Require: Number of unbroken rays n u . 

{Algorithm for reconstructing f(x,y) in £1 in the presence of obstacle Oo} 
initialize an empty list L that will contain all broken and unbroken rays 
for p = 1 — >■ n& do 

generate a random broken ray r from input data R, T, H, 0,q and Oi {a broken ray is a 
unique triple ((x r , y r ), (x t , yt), (x^, Vh)) that does not intersect the obstacle except at the 
reflection point (xh,yh)} 
add r to L 
end for 

for p = 1 — >■ n u do 

generate a random unbroken ray r from input data R, T, H, £Iq and £l\ {an unbroken ray 
is a unique pair ((x r ,y r ), (x t ,yt)) that does not intersect the obstacle} 
add r to L 
end for 

randomize L as a preprocessing step before starting the Kaczmarz method 
for all rays r in L do 

Compute travel time p r for ray r. In numerical simulations, this is obtained via nimerical 
integration in fio along r. Store p r in list P. 
end for 

Reconstruct f(x,y) by the Kaczmarz method [11] for the linear system associated with rays 
r in L and corresponding travel times p r in P. 

As input to the above algorithm give ra& and n u to be equal or approximately equal to 
the maximum number of broken and unbroken rays for the finite input sets R, T and H, or 
alternatively, during ray generation, generate new rays until all rays with end points in the 
input sets are generated. In addition, again in order to approximate the requirements of the 
theory of broken ray tomography for inclusion of all broken and unbroken rays, provide as input 
as many transmitter and receiver points on the observation boundary as possible. Ideally, the 
transmitter and receiver points are uniformly distributed on the boundary and sufficiently close 
to each other so that each set approximates the set of all points on the observation boundary. 
The above algorithm uses Lambertian reflection by considering broken rays as random triples 
((x r ,y r ), (x t ,yt), (xh,Uh)) and uses all rays in the domain. 



3. Experimental Results 

In order to show the effectiveness of the numerical solution of the broken ray tomography 
problem with Lambertian reflection, a Java implementation of the above algorithm compares 



the reconstructed values of f(x,y) in f^o with the known values of f(x,y) for the same test 
function f(x,y) = K\J (x — xn) 2 + {y — yo) 2 , where (xo,yo) is the center of the computation 
grid, and test environment as in [10]. The following table compares the reconstruction error for 
classical tomography without reflection and broken ray tomography with Lambertian reflection. 



Experiment 


ART Error 


ART Iterations 


BRTL Error BRTL Iterations 


1 


1.269592e-004 


54293 


1.272610e-005 


32635 


2 


1.300814e-004 


54204 


5.272083e-006 


40557 


3 


1.478464e-004 


45923 


7.530267e-006 


37251 


4 


1.642026e-004 


42480 


1.110909e-005 


33104 


5 


1.927561e-004 


23037 


1.163345e-005 


33424 


6 


1.985439e-004 


33705 


1.896151e-005 


28522 


7 


8.991362e-005 


88190 


1.401511e-005 


30061 


8 


1.641201e-004 


40717 


1.955676e-005 


29687 


9 


1.089771e-004 


64932 


7.102756e-006 


40126 


10 


1.934379e-004 


32899 


2.149935e-005 


28723 


Average 


1.516838e-004 


48038 


1.294065e-005 


33409 



Table 1. Error and number of iterations for broken ray tomography with Lambertian 
reflection(BRTL) at the boundary of the reflecting obstacle for a fixed number of 126050 rays. 
The average error for BRTL is 1.294065e-005 and the average number of iterations for finding a 
solution is 33409, and are shown in the right two columns of the table. The results for classical 
tomography with the ART method with the same number of rays are shown in the left two 
columns of the table. The average error for ART is 1.516838e-004 and the average number of 
iterations is 48038. 

The results show that the reconstruction error of the new numerical solution of the broken ray 
tomography problem using Lambertian reflection is approximately three times smaller compared 
to an average reconstruction error using specular reflection of 3.874848e-005 by the algorithm 
from [10] for the same test function, domain and obstacle and an order of magnitude smaller 
than the reconstruction error for classical tomography in the presence of a reflecting obstacle. 
The new method also appears to be faster compared to classical tomography with ART. The 
order of the rays in both the ART and BRTL tests in this work is randomized before applying the 
Kaczmarz method therefore the speed of convergence difference is due to the use of reflection. 
We will report further results on the speed of convergence. 
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